SLAC-PUB-12509 
May, 2007 



Studies of Cosmic Rays with GeV Gamma Rays * 



Hiroyasu Tajima,^ Tuneyoshi Kamae/ Stefano Finazzi,^ Johann Cohen- Tanugi^ and James Chiang^'^ 

(1) Stanford Linear Accelerator Center, 2575 Sand Hill Road, Menlo Park, CA 94025, USA 
(2) Scuola Normale Superiore, 1-56100 Pisa, Italy 
(3) Center for Research and Exploration in Space Science and Technology (CRESST), 
University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA 



Abstract 

We describe the role of GeV gamma-ray observations with GLAST-LAT (Gamma-ray Large Area 
Space Telescope - Large Area Telescope) in identifying interaction sites of cosmic-ray proton (or 
hadrons) with interstellar medium (ISM). We expect to detect gamma rays from neutral pion decays 
in high-density ISM regions in the Galaxy, Large Magellanic Cloud, and other satellite galaxies. These 
gamma-ray sources have been detected already with EGRET (Energetic Gamma Ray Experiment 
Telescope) as extended sources (eg. LMC and Orion clouds) and GLAST-LAT will detect many more 
with a higher spatial resolution and in a wider spectral range. 

We have developed a novel image restoration technique based on the Richardson- Lucy algorithm 
optimized for GLAST-LAT observation of extended sources. Our algorithm calculates PSF (point 
spread function) for each event. This step is very important for GLAST-LAT and EGRET image 
analysis since PSF varies more than one order of magnitude from one gamma ray to another depending 
on its energy as well as its impact point and angle in the instrument. The GLAST-LAT and EGRET 
image analysis has to cope with Poisson fluctuation due to low number of detected photons for 
most sources. Our technique incorporates wavelet filtering to minimize effects due to the fluctuation. 
Preliminary studies on some EGRET sources are presented, which shows potential of this novel image 
restoration technique for the identification and characterisation of extended gamma-ray sources. 
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1. INTRODUCTION 



The origin of the cosmic rays has been a great mystery since their discovery by Victor Hess in 1912. Supernova 
remnants (SNRs) are considered to be the best candidates for the proton (hadron) acceleration up to the so-cahed 
"knee" (3 x 10^^ eV) in the cosmic-ray spectrum. Supports for this hypothesis are mostly circumstantial and 
theoretical rather than observational. For example, the kinetic energy released by supernova explosion has been 
known to be comparable with the total energy in the Galactic cosmic ray.^ Simple diffusive shock acceleration 
modeling of a typical SNR has given the maximum energy of accelerated particles to reach the " knee" (see for 
example a review article^). Synchrotron emission in the X-ray band observed in SN1006 by ASCA^ was the first 
strong indirect evidence for the existence of ^100 TeV electrons in an SNR. Observations of TeV gamma rays 
from another SNR, RX J1713.7-3946, by CANGAROO^ and by H.E.S.S.^ were the first confirmed direct evidence 
for particle acceleration up to 100 TeV in a SNR. Furthermore, H.E.S.S. provided the first gamma-ray image of 
RX J1713. 7-3946, which is qualitatively consistent with that obtained in the X-ray band. This morphorogical 
coincidence has been interpreted that the accelerated particles are interacting with ISM or cosmic microwave 
backgrounds (CMB) near or in the shock region. However, the TeV spectrum observed by H.E.S.S. alone does 
not favor decisively a dominant gamma-ray emission mechanism (Compton up-scattering or 7r^-decays following 
proton-proton interactions). Hence, the proof of hadronic acceleration in SNRs is not yet conclusive and wait 
for new data from GLAST-LAT in the GeV energy range. 

There is no consensus for acceleration sites of cosmic rays above the "knee" of the cosmic-ray spectrum. 
Several astronomical objects such as gamma-ray bursts (GRBs), active galaxies or merging galaxy clusters have 
been put forward as possible candidates. Observation of UHECRs (Ultra High-Energy Cosmic Rays) and UHE- 
neutrinos associated with known astronomical sources will provide convincing evidence for UHECR acceleration 
in these sources. Unfortunately the phase space for UHECRs to reach a detector is rather small because charged 
particles are deflected by the extragalactic and interstellar magnetic field below 10^^ eV. It is also known that 
hadronic cosmic rays cannot travel more than ~20 Mpc above ~ 10^^*^ eV due to production of the A resonance 
by interacting with CMB (known as GZK suppression^). Neutrinos are free from both constraints and allow us 
to probe distant UHECR acceleration sites if there is enough target material at the sites. However, neutrinos 
are difficult to detect due to their small cross section. In the case of gamma-ray observations, mult i- wavelength 
analysis is critical to resolve between the leptonic and hadronic model of gamma-ray emission. Neutrino detection 
from the same gamma-ray sources wil provide definitive evidence for acceleration of protons in such sources. 

2. COSMIC RAYS FROM SUPERNOVA REMNANTS 

As described above, SNRs are the prime candidates for the origin of Galactic cosmic rays. Recent observations of 
TeV gamma rays from RX J1713. 7-3946 and RX J0852.0-4622, by H.E.S.S.^' ^ present undisputed evidence that 
electrons and/or protons are accelerated to at least 10"*^^ eV, which is very close to the "knee" energies. Fig. 1 (a) 
shows the latest TeV gamma-ray image of a shell-type SNR, RX J1713. 7-3946, which is the first SNR imaged 
in TeV gamma rays. Morphological similarity between X-ray^ and TeV gamma ray^ as shown in Fig. 1 (a) 
may suggest that Compton up-scattering of CMB photons to TeV by high energy electrons and positrons that 
produce the synchrotron radiation observed by ASCA. The azimuthal profile of gamma-ray intensity observed 
by H.E.S.S. and CO intensity observed by NANTEN are not as weU correlated as between H.E.S.S. and ASCA. 
This may present some difficulties for the proton interaction model. The spectrum of the H.E.S.S. observation 
does not necessarily agree with the prediction by the leptonic model. Furthermore, Hiraga et al}^ pointed 
out that the radial profile of the X-ray mission cannot be explained by smooth spherical distributions expected 
from leptonic models. In addition, they observed possible evidence of the positive correlation between the 
X-ray brightness and the absorption column density, which may indicate that molecular clouds are providing 
unaccelerated electrons (and protons) to the shock, resulting in higher electron population thus higher X-ray 
emissivity. Using this scenario, the morphological similarity between TeV gamma ray and X-ray images can be 
explained as a similarity in proton and electron injection rate in the dense matter region, making a hadronic 
model a viable alternative to leptonic models. More definitive evidence is required to conclusively determine 
the dominant gamma ray emission mechanism at RX J1713. 7-3946. Fig. 1 (b) shows the H.E.S.S. measurement 
result of RX Jl 713. 7-3946 along with simulated GLAST-LAT measurements for hadronic and leptonic cases. 
It demonstrates that GeV gamma-ray observations by the GLAST-LAT can provide conclusive evidence on this 
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matter since the spectra predicted by hadronic and leptonic models are clearly separated in the GLAST-LAT 
range (-100 MeV to 300 GeV). 




3. COSMIC-RAY INTERACTIONS WITH THE GALACTIC INTERSTELLAR 

MEDIUM 

Diffuse gamma-ray emissions from the Galactic plane provide rich information on Galactic cosmic rays, cosmic-ray 
propagation, Galactic interstellar medium (ISM), and cosmic-ray interactions with ISM and low-energy photons. 
Hunter et al}^ fitted the EGRET gamma-ray spectrum in the Galactic center region with a model gamma-ray 
spectral shape obtained on a model proton-proton pion production cross-section and a model for the spatial 
distribution of pion gamma-rays projected to the Galactic coordinate. The fit left the absolute normalization 
undetermined because matter density nor cosmic-ray intensity in the Galactic center region were not known. The 
spectral shape of the model used in the publication reproduced the observed spectrum poorly when normalized 
at the peak in the vF{v): the model left possible 50% "excess" in the observed flux above 1 GeV. This possible 
excess gamma-ray flux has led to many speculations including claims for dark matter detection. Mori^^ fitted the 



(a) (b) 




Figure 2. (a) The diffuse gamma-ray differential spectrum, multiplied by for 330° < £ < 30°, |6| < 6°.-'^^ Data points are from 
EGRET observations. Model calculations are superimposed by curves, (b) EGRET diffuse gamma-ray image in the Galactic center 
region (top panel) and CO survey map in the same region-*^^ (bottom pannel). 
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observed spectrum with a simulation code (Pythia) developed for particel physics experiments. He noted that 
this excess can be explained by a harder cosmic-ray spectrum, 1.48 x ^~2.45 cm~^s~^sr~^GeV~^. He suggested 
that the cosmic-ray spectrum in the Galactic center region could be harder than that in the solar neighborhood, 
which was used for the above model calculations. Kamae et al^^ corrected the old nuclear interaction models 
and demonstrated that 20-50% more pions are produced compared with old models used by Hunter et al by 
correctly accounting for the diffraction dissociation, scaling violation and rising inelastic cross-section in the 
nuclear interactions as shown in "Model A with LIS (Local Interstellar Spectrum)" of Fig. 2 (a). However, the 
remaining excess still requires a harder cosmic-ray spectrum as shown in "Model A with Trial4GR" of Fig. 2 (a) 
where the power law index is set to —2.5 instead of —2.7 above 20 GeV. 

Obviously, better understanding of the cosmic-ray spectrum in the Galactic center region is crucial to solve 
this problem. In addition, contamination by unresolved gamma-ray point sources may complicate this problem. 
Systematic studies of high-density regions well surveyed by CO observations (both inside and outside of the 
Galactic plane) and associated gamma-ray emissions are extremely useful. For some cases the CO observation 
can provide 3D information of the ISM, which can be used in modeling acceleration and propagation of protons 
and electrons at selected supernova remnants. Fig. 2 (b) shows the EGRET diffuse gamma-ray image in the 
Galactic center region (top panel) and the CO map in the same region^^ (bottom pannel). For dense molecular 
clouds, CS survey^'' and other CO transitions can supplement the CO data. Finer and more complete surveys 
(eg. Ref. 18) are becoming available in this region. 

4. COSMIC RAYS FROM EXTRA-GALACTIC SOURCES 

Gamma-ray bursts (GRBs) are among the most energetic explosions in the Universe, originating possibly in core- 
collapse of massive stars or binary mergers of neutron stars and/or black-holes. There are theoretical speculations 
that cosmic rays can be accelerated beyond the "knee" in GRBs. We expect that the prompt gamma-ray emissions 
come from compact regions because of the observed variability (down to milliseconds). The gamma-ray power 
spectrum of the GRBs is believed to follow broken power law with the peak energies at around 0.1-1 MeV. This 
spectral feature is interpreted as due to synchrotron radiation or Compton up-scattering in highly relativistic 
internal shock regions. However, the delayed high energy gamma-ray emission {E > 10 MeV) observed in the 
GRB941017^^ as shown in Fig. 3 (a) presents a serious challenge to the standard interpretation. Soft gamma-ray 
emission decayed in 100 s while the high energy component remained stable. Such difference in time structure is 
difficult to explain with simple electron models. The photo-hadronic (^7) interactions may produce GeV-TeV 
gamma-ray emissions due to secondary electrons from hadronic and electromagnetic cascades. GLAST-LAT 
will observe many more GRBs like GRB941017 and clarify the origin of high-energy gamma rays. Such data 



(a) (b) 




Figure 3. (a) Time dependent energy fluxes for GRB941017.19 Crosses and circles correspond to BATSE-LAD and EGRET-TASC 
data respectively, (b) Model prediction of gamma-ray spectrum for a merging galaxy cluster, Coma.^-*^ Expected sensitivity limit 
of GLAST-LAT and lACTs are also displayed. 
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will shed new light into understanding the GRB mechanism. In particular, finding out the peak energy and the 
power of the high energy component will provide strong contraints on GRB gamma-ray emission models. 

Merging galaxy clusters are also considered to be one of candidates for the origin of UHECRs because the 
Lamor radius can reach Mpc. UHECRs (protons) accelerated in such acceleration sites will interact with CMB 
to producing UHE electrons and positrons over cosmological time scale. These secondary leptons can produce 
hard X-rays via synchrotron radiation and gamma rays via Compton up-scattering of CMB photons. Multi- 
wavelength observation of merging galaxy clusters in X-ray and gamma-ray bands can constrain such models as 
demonstrated by Fig. 3 (b).^^ It shows a model prediction of X-ray and gamma-ray spectra with upper limits 
given by existing X-ray observations for the Coma galaxy cluster. Observation of the Coma with GLAST-LAT 
(and Imaging Atomospheric Cerenkov Telescopes in the northern hemisphere such as MAGIC or VERITAS) will 
either confirm or put a stringent constrains to this scenario. 



5. EXTENDED SOURCE ANALYSIS FOR GLAST 

The sources discussed above are most likely spatially extended. Modest PSF (point-spread function) and photon 
statistics of the GLAST-LAT compared with X-ray instruments present challenge to studies on these extended 
sources. Morphological feature of TeV sources has played a critical role in identifying all SNRs and PWNe 
(Pulsar Wind Nebulae) in H.E.S.S. analyses. Several TeV sources remain unidentified probably due to limited 
spatial resolution in the TeV band. The PSF of GLAST-LAT is often comparable to size of extended sources 
and morphology of sources may not be evident. In addition, the presence of the Galactic diffuse gamma-ray 
background makes it even more difficult to resolve extended sources. The fact that roughly half of 271 gamma-ray 
sources^^ detected by EGRET remain unidentified underscores the challenge GLAST-LAT will face. We note 
that many of the unidentified EGRET sources are located at low Galactic latitude and are considered to be of 
Galactic origin. Furthermore, extended sources will be contaminated with contribution from nearby sources. 

The iterative image restoration technique originally suggested by Richardson and Lucy^^'^^ has potential 
to improve the quality of GLAST-LAT and EGRET images. Hook and Lucy also developed the dual channel 
method^^ to account for the effect of point sources in the restoration process, which is effective in removing 
point-source contibution from a diffuse image. However, this technique tends to suffer from amplification of the 
Poisson noise as the number of iteration is increased. Efforts have been made to define objective criteria to 
stop the iteration before such effects become prominent. Charalabides et al?^ applied this technique to the 
EGRET images using a wavelet filtering to suppress the effect of Poisson noise. However, the results suffered 
from rapidly- varying and energy-dependent PSF. We have developed a novel image restoration technique that 
incorporates PSF calculated for each event to address this problem. This technique will be critical for the analysis 
of extended sources in the GLAST-LAT since its PSF varies by two orders of magnitude over the wide energy 
band, 0.1-300 GeV. 

5.1. Image Restoration Technique 

The observed image ^{x) is a convolution of the true image '^(^) and the instrument response function P(x|^), 
where P{x\£^) is the probability that the photon is observed at x when the true position is ^. Based on a Bayesian 
approach, the original image '0(C) can be obtained iteratively: 

It is proven mathematically that yields larger likelihood than '0^(0 at each iteration. We use wavelet 

denoising to suppress the noise in the residual between the observed and expected image in each iteration. 
The residual, p^{x)^ is defined as p^{x) = — where (jf{x) is the convolved image: 

clf{x) = J P{x\OriOdC (2) 
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The residual is decomposed by the "a trous" wavelet algorithm into J components: 

J 

(x) = cj {x) + ^ Wj {x) , (3) 

where cj{x) is the last smoothed image and Wj denotes the wavelet component at scale j. Denoising is performed 
by requiring the values to be above noise in each component: the filtered residual is written as 

J 

r{x) = cj{x) -^^Mj{x)wj{x) , (4) 

where Mj{x) is if Wj{x) is consistent with noise and 1 otherwise (hard-thresholding). We currently use the 
99% Poisson probability as the consistency threshold. The filtered residual is incorporated in the deconvolution 
process by replacing the observed image 0(x) = + with the filtered image (x) = + (t>'^{x)^ 



■'«'=/ 



When we have known point sources, the probability density function, ^^(^), can be decomposed into two 
components: = ^l^ps{^,) + V^dfIO^ where ?/^p^(^) represents a point source component (mostly zero apart 

from ^-functions at the positions of the point sources) while '0df(O represents a diffuse component. Wavelet 
filtering is applied only to the diffuse component. This dual channel method enables us to incorporate sharp 
point sources while keeping the smooth diffuse component. This effectively works as point source subtraction 
without negative values in '^^^(O (straight point source subtraction often causes negative values due to noise). 

This image restoration technique can be generalized for the case where the PSF is assigned for each event as: 

where Xk is the observed position of the kth event and Pk{xk\(,) is the probability to observe the event at Xk 
when the true position is ^. In this case, the convolved image is defined as 

cp^ix) = I UixlOriOdC, (7) 

where n(x|C) is the weighted average of the instrument response at a given position x and defined as: 

n(.|c) = (8) 

The resulting residual, p^{x), can be filtered with the same wavelet denoising method as above to obtain p^{x). 
The filtered image cj) (x) =]f{x) + ^''(x) can be incorporated in the iteration process as: 

These generalized formula described here are completely equivalent to the non- generalized version if Pk{xk\£,) is 
common for all events. 
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5.2. Application to EGRET Image Analysis 

In this section, we demonstrate our image restoration technique by analyzing three fields observed by EGRET. 
The PSF is calculated for each event as a convolution of PSFs in different energy bins, according to the energy 
dispersion of the given energy measurement. Only photons with energy greater than 0.5 GeV are used in this 
analysis in order to minimize the processing time without sacrificing the image quality. (PSF in this energy region 
varies from ~ 2° to 0.5°.) Lower-energy photons do not contribute to improving the quality of the restored image. 
It does not preclude the use of lower-energy photons in the future analysis, in particular for spectral analysis, 
since it does not degrade the image quality either since such photons are weighted appropriately via a large PSF. 
The results presented in this paper are still preliminary and should not be used for any quantitative analysis. 

5.2.1. Large Magellanic Cloud 

The LMC (Large Magellanic Cloud) is one of few extended sources resolved by EGRET. Fig. 4 shows a smoothed 
photon count map of the LMC in the left panel, a restored image of the LMC in the middle panel, and the IRAS 
LMC image overlaid with the contour of the restored EGRET image in the right panel. The restored EGRET 
image is apparently sharper with better contrast than the original count map while avoiding enhancement of 
spurious structures (sometimes suppressing them) in the background (The color scale is common for both images). 
The bright spot of the restored image is well correlated with a dense region indicated by the IRAS image, which 
traces dusty interstellar clouds and associated star-forming regions. 




Figure 4. Smoothed count map (left), deconvolved image (middle) of the LMC by EGRET, and the IRAS LMC image overlaid 
with contour of the deconvolved EGRET image (right). 

5.2.2. 3EG J1234-1318 

We also analyzed several unidentified EGRET sources at high Galactic latitude, where the Galactic diffuse 
background is very low. We found one source that is inconsistent with one point source as shown in Fig. 5 The 
left panel shows the smoothed count map of the 3EG J1234-1318 while the right panel shows the restored image. 
A white cross and a contour in each image indicate the source location and its 90% confidence level contour in 
the original 3EG catalog (The 3EG catalog is produced using photons with E > 0.1 GeV). The restored image 
has much better contrast and clearly shows an elongated structure indicating the presence of either an extended 
source, or two point sources separated by ~1.5°. Since the Galactic latitude of this EGRET signal is ~50° 
it is very likely to be two extra- Galactic point sources. If such is the case, the preliminary location of a new 
EGRET point source would be (RA, DEC) = (187.9°, -14.1°) for the brighter one and (188.9°, -13.0°) for the 
fainter one instead of the original position, (188.5°, —13.3°). Another possibility is that this is due to a nearby 
molecular cloud interacting with cosmic rays, which could be extremely useful to obtain cosmic-ray fiux (both 
spectrum and intensity) in that region. This result demonstrates that our image restoration technique can be 
quite effective in identification of extended sources and source confusion. 
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Figure 5. Smoothed count map (left) and deconvolved image (right) of the source 3EG J 1234-1318 by EGRET. 

5.2.3. Galactic Plane 

Although the Milky Way can be resolved by the EGRET due to its large size, smaller scale structure in the 
Galactic plane is not easy to resolve due to the large EGRET PSF and contamination by the point sources. Our 
image restoration technique can reduce these effects as shown in Fig. 6 The top left panel shows the smoothed 
count map, the top right panel shows the simple deconvolved image (no dual channel method), the bottom left 
panel shows the deconvolved image with known point sources included by the dual channel method, and the 
bottom right panel shows the deconvolved image of the only diffuse channel from the previous image, which is 
equivalent to the point-source removal. The final EGRET image of the Galactic plane has smooth background 
and less structures associated with "point sources" that may or may not be true point sources. 




o do rio\ 



10° 5 



335° 350° 



do rio\ 



10° 5 



335° 350° 



Figure 6. Smoothed count map (top left), simple deconvolved image (top right), deconvolved image with known point sources 
included by the dual channel method (bottom left), and deconvolved image of the only diffuse channel (bottom right) of the Galactic 
plane around the Galactic center by EGRET. 
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6. CONCLUSIONS 



GeV gamma-ray observations by GLAST-LAT will shed light on the origin, propagation and interaction with 
ISM of cosmic rays. Specifically, the GLAST-LAT will conclusively distinguish leptonic and hadronic models 
of gamma-ray emissions in SNRs RX J1713. 7-3946 and RX J0852.0-4622. The GLAST-LAT may give clues on 
understanding the cosmic-ray acceleration in extra- Galactic sources such as GRBs and clusters of galaxies. 

In order to facilitate the analysis of extended sources with limited PSF expected for the GLAST-LAT, 
generalization of a widely-used image restoration technique is introduced to incorporate the strongly energy- 
dependent PSF in the EGRET and GLAST-LAT although it is not trivial since the wavelet filtering is performed 
on pixelated images. We demonstrated the effectiveness of our technique using three fields observed by EGRET. 
Our technique is confirmed to be useful to obtain sharper image of extended sources, to identify source confusion, 
and to cleanly remove point sources from the image of extended sources. 

This work was supported in part by the U.S. Department of Energy under Grant DE-AC02-76SF00515. 
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